A Markov chain is a stochastic process \(\{X_t\}\) indexed by time \(t\geq 0\).

• We call the set of all values that \(\{X_t\}\) can assume the state space of the Markov Chain.
• \(\{X-t\}\) is called a Markov Chain if it fulfills the Markov property that \[ P(X_{t+1}=j\mid X_0=i_0,\dots,X_t=i_t)=P(X_{t+1}=j\mid X_t=i_t), \] i.e., that the conditional distribution of the next state depends only on the current state.

Markov Chain Monte Carlo (MCMC) methods make use of Markov chains to sample from a target distribution.

In non-physical applications,

The system is described by a so-called Hamiltonian equation, \(H(q_p)\). Its partial derivatives, the equations of motion, determine how \(\boldsymbol{q}\) and \(\boldsymbol{p}\) change over time:

\[ \begin{aligned} \frac{\mathrm{d}q_i}{\mathrm{d}t}&=\frac{\partial H}{\partial p_i},\\ \frac{\mathrm{d}p_i}{\mathrm{d}t}&=-\frac{\partial H}{\partial q_i},\\ \end{aligned} \]

for \(i = 1, \dots, d\), and \(d\) being the length of the vectors \(\boldsymbol{q}\) and \(\boldsymbol{p}\).

For Hamiltonian Monte Carlo, we usually use the following kind of Hamiltonian functions:

\[ H(\boldsymbol{q},\boldsymbol{p}) = \textcolor{var(--secondary-color)}{U(\boldsymbol{q})} + \textcolor{var(--tertiary-color)}{K(\boldsymbol{p})}, \]

where the potential energy \(U(\boldsymbol{q})\) is defined as minus the log probability density of the distribution for \(\boldsymbol{q}\) that we wish to sample, plus a constant;

and the kinetic energy \(K(\boldsymbol{q})\) is given by

\[ \textcolor{var(--tertiary-color)}{K(\boldsymbol{p})} = \boldsymbol{p}'\boldsymbol{M}^{-1}\boldsymbol{p}/2, \]

where \(\boldsymbol{M}\) is a symmetric, p.s.d., and often diagonal mass matrix.

Hamiltonian dynamics fulfill a set of properties that make them suitable for use in MCMC updating:

• Reversibility. The mapping from the state at time \(t\) to the state at time \(t+s\) is one-to-one and hence has an inverse. This means that Markov chain transitions are reversible.
• Conservation of the Hamiltonian. With dynamics as specified, the Hamiltonian \(H\) itself is kept invariant. For Metropolis updates, the acceptance probability is one if \(H\) is kept invariant (only approximatively achievable).
• Volume preservation. Hamiltonian dynamics preserves volume in \((\boldsymbol{q},\boldsymbol{p})\) space, meaning that we do npt need to account for a change in volume in the acceptance probability for Metropolis updates.

Reversibility and preservation of volume can be maintained even when the Hamiltonian is approximated.

For implementation, we need to discretize the Hamiltonian equations using a small step size \(\varepsilon\): time is then discrete with \(t = 0, \varepsilon, 2\varepsilon, 3\varepsilon, \dots\)

The simplest way to approximate the solution is Euler’s method:

We can obtain better results by slightly modifying Euler’s method:

For even better results, we can use the Leapfrog method:

Using Hamiltonian dynamics to sample from a distribution requires translating the density to a potential energy function and introducing momentum variables. We can use the concept of a canonical distribution from statistical mechanics. Given an energy function \(E(x)\) for a state \(x\), the canonical distribution over states has density \(P(x)=(1/Z)\mathrm{exp}(-E(x)/T)\), where \(T\) is temperature and \(Z\) is a normalizing constant. Using the Hamiltonian \(H(\boldsymbol{q},\boldsymbol{p})=U(\boldsymbol{q})+K(\boldsymbol{p})\) as an energy function, we get

\[ P(\boldsymbol{q},\boldsymbol{p}) = \frac{1}{Z}\mathrm{exp}(-U(\boldsymbol{q})/T)\mathrm{exp}(-K(\boldsymbol{p})/T). \]

We can see that \(\boldsymbol{q}\) and \(\boldsymbol{p}\) are independent, and both have canonical distributions. The former will be used to represent our variables of interest, and the latter for momentum.

Each iteration has two steps:

1. In the first step, new values for the momentum variables are drawn from their Gaussian distribution.
2. In the second step, a Metropolis update is performed using Hamiltonian dynamics as explained before. We get a proposed new state \((\boldsymbol{q}^*,\boldsymbol{p}^*)\). This state is accepted as the next state of the Markov chain with probability \(\mathrm{min}\Big(1,\quad \mathrm{exp}(-H(\boldsymbol{q}^*,\boldsymbol{p}^*)+H(\boldsymbol{q},\boldsymbol{p}))\Big)\quad =\quad \mathrm{min}\Big(1,\quad \mathrm{exp} \big(-U(\boldsymbol{q}^*)+U(\boldsymbol{q})-K(\boldsymbol{p}^*)+K(\boldsymbol{p})\big)\Big).\) If the proposed state is not accepted, the next state is set equal to the current state.

Only in the very first step of the chain does the probability density for \((\boldsymbol{q},\boldsymbol{p})\) change from one step to the next. Also, the algorithm leaves the canonical distribution invariant.

We can use HMC to sample only from continuous distributions on \(\mathbb{R}^d\) for which the density function can be evaluated and the partial derivatives of its log can be computed.

It allows us to sample more efficiently from such distributions than simpler methods such as random-walk Metropolis.

However, in order to actually attain the benefits of HMC, we have to properly tune the step size \(\varepsilon\) and the trajectory length \(L\), i.e., the number of leapfrog steps. Tuning an HMC algorithm is also more difficult than tuning a simple Metropolis method.